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Abstract 

We study the dynamics of coherent waves in nonhnear honeycomb lattices and show that nonlin- 
earity breaks down the Dirac dynamics. As an example, we demonstrate that even a weak nonlin- 
earity has major qualitative effects one of the hallmarks of honeycomb lattices: conical diffraction. 
Under linear conditions, a circular input wave-packet associated with the Dirac point evolves into 
a ring, but even a weak nonlinearity alters the evolution such that the emerging beam possesses 
triangular symmetry, and populates Bloch modes outside of the Dirac cone region. Our results are 
presented in the context of optics, but we propose a scheme to observe equivalent phenomena in 
Bose-Einstein condensates. 
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The past decade has witnessed considerable interest in honeycomb lattices in many fields, 
ranging from carbon nanotubes l| and graphene in condensed matter, cold atoms in 

loneycomb optical lattices js-jsl, and electromagnetic waves in honeycomb photonic crystals 
9l-|l2j and photonic lattices (waveguide arrays) 13l-ll5|. The unique features of the honey- 
comb lattice result from its unusual band structure, displaying two intersecting bands. The 
vicinity of the intersection points is conical, and excitations residing at the close vicinity 
of the intersection points obey the massless Dirac equation: a relativistic wave equation 
for massless spin-half particles. As a consequence of this unusual linear dispersion and the 
additional degrees of freedom (two sites in each unit cell, referred as pseudo-spin), a vari- 
ety of exciting phenomena have been obtained: room temperature quantum Hall effect in 
graphene 16|], conical diffraction in hoiieycomb photonic lattices 



sistance which implies anti- localization 17H19|. Klein tunneling 
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ISL negative magnetore- 
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2211. and Zitterbewe- 



gung |2J] to name a few. The dynamics of Dirac-like excitations in honeycomb lattices 
has been well studied when the propagation equation is linear (interaction-free). However, 
the dynamics can also be nonlinear, as happens in the optical domain due to light-matter 
interactions, and for atomic Bose-Einstein condensates (BEC) due to pair-wise scattering. 
In either of these, the nonlinear dynamics has received little attention in the context of hon- 
eycomb lattices. The first study of nonlinear dynamics in honeycomb lattices was conducted 
in 13[, demonstrating gap solitons, which had no overlap with Bloch modes residing at the 
vicinity to the Dirac points. Later studies of nonlinear dynamics in honeycomb lattices were 
a generalization of Dirac approximation 25|, |26| , where a nonlinear version of the massless 



Dirac equation have been studied. 

Here, we study the nonlinear dynamics of waves in honeycomb lattices, where initial 
wave-packets are comprised of Bloch waves from the vicinity of the Dirac points. Such wave- 
packets are very well described by the massless Dirac equation. However, even the presence 
of a fairly weak nonlinearity drives the waves away from the vicinity of the Dirac points, 
hence the nonlinearity breaks the Dirac dynamics associated with honeycomb lattices. We 
demonstrate this dramatic change in the dynamics by investigating conical diffraction: an 
initial wave-packet with circular symmetry that evolves into two concentric rings (separated 
by a dark ring). Surprisingly, under slightly nonlinear conditions, the same (circular) input 
beam evolves into a triangular-ring beam. Interestingly, we find that the resulting triangle 
is rotated by vr when the nonlinearity changes sign, i.e., a given input beam subjected to a 
focusing nonlinearity (attractive interactions for BEC) evolves into a triangle that is rotated 
by TT with respect to the triangle obtained when the same input beam is subjected to a de- 
focusing (repulsive) nonlinearity. Hence, the presence of nonlinearity (interactions) changes 
the wave dynamics in honeycomb lattices dramatically, introducing important changes to the 
unique phenomena associated with the effective Dirac equation, such as conical diffraction 
and Klein tunneling. Moreover, the nonlinear breakdown the Dirac dynamics - in itself gives 
rise to interesting new effects. Our results are presented in the context of optics, however, 
we propose a scheme to observe equivalent phenomena in Bose-Einstein condensates. 
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The paraxial propagation of a monochromatic field envelope ip inside a photonic lattice 
with Kerr non-linearity is described by 

where 6n{x, y) is the modulation in the refractive index defining the lattice as shown in Fig. 
[1] (a), k is the wave-number in the medium, rig the background refractive index, and n2 
is the Kerr coefficient. The sign of n2 determines the type of nonlinearity, where n2 > 
corresponds to a focusing nonlinearity (attractive interactions in the context of BEC). The 
term k6n/n0 is referred to as the optical potential. 

At low intensities (n2|^/'P <^ max{6n{x,y)}), the nonlinear term is negligible and the 
propagation can be described by linear methods. In the absence of the nonlinear term, 
Eq.([T]) has solutions of the form ip{x,y, z) = U{x,y) exp{if3z) where U{x,y) is a solution of 
the eigenvalue equation 

HoUix,y) = pUix,y), where Ho = -^Vl - —6nix,y), (2) 

and the eigenvalue, (3, is the propagation constant (analogous to the energy in quantum 
mechanics with opposite sign, i.e., (3 ^ —E). Since 6n{x,y) is a periodic function in x and 
y, the eigenfunctions of Hq can be chosen to be Bloch waves 

Bq,m{r) = Uq^rn{r) cxp {iq ■ r), (3) 

where Uq^n{x, y) has the same periodicity as the potential, q = (g^., qy) is the lattice momen- 
tum, r = [x, y) is the transverse coordinate, and m is the band index. We solve numerically 
Eq.([2]) in the first Brillion zone (BZ) (Fig. [1] (b)) and obtain the band structure (3{q) (Fig. 
[It^c)). Note the intersection of the bands at the corners of the first BZ. According to Fig. 
[Dj^a), it appears as if there are six intersection points; however, these six points are in fact 
two triplets of two non-equivalent points denoted by K±. (Fig. [1] (b)). That is, there are 
only two non-equivalent intersection points, since the corners are connected by reciprocal 
lattice vectors (Fig. [11(b)). The close vicinity of the intersection points is conical (Fig. [11(d)). 
In the tight binding approximation, the effective Hamiltonian describing the Bloch modes 



from the vicinity of these points is the Dirac Hamiltonian for massless particles 27|. Hence, 
these points are often referred to as Dirac points, and the band structure in the vicinity of 
these points forms Dirac cones. We emphasize that Bloch modes that reside outside the 
Dirac cone cannot be described by the Dirac equation, since the propagation constant (3{q) 
is no longer linear in q, and moreover it does not possess circular symmetry but rather a 
three-fold symmetry (as will be discussed later). 

Our main interest in the article is the effect of nonlinearity on the extremely unique 
diffraction pattern, named conical diffraction, obtained for an input beam comprised of 
Bloch modes residing in the Dirac cone. Conical diffraction refers to the fact that at a large 
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enough propagation distance, the circular (envelope) input beam evolves into two bright 
ring s of constant width separated by a dark ring, and a vanishing intensity at the origin 
A characteristic conical diffraction pattern is shown in Fig. [2](b), where the first 
ring is separated from the second ring that is just starting to form (for two distinct rings 
see Experimentally, the input beam can be constructed by means of a spatial light 

modulator (SLM), or by three plane waves with wave- vectors corresponding to the three 
equivalent corners (say Kj^) of the BZ 141 ] . 

At high intensities, the nonlinear term comes into play. The input wave-packet has an 
envelope that possesses circular symmetry. Hence, the optical potential induced through 
the nonlinear term has circular symmetry as well. The initial excitation is around the Dirac 
points, where the Dirac cone approximation is excellent. Thus, one could expect that, even 
though some qualitative differences would emerge, the resulting diffraction pattern will still 
have circular symmetry. This educated guess is supported by analysis of the "nonlinear Dirac 
equation" (Dirac equation with a nonlinear term) 25|. We test this idea by solving Eq.([T]) 
numerically using the split-step Fourier algorithm in the presence of a focusing nonlinearity, 
where the input beam is comprised of Bloch modes from the vicinity of the point K^. We 
construct the input beam using k ■ p approximation 

^(r, 0) = VJ7[BK^^^{r) + zfiK+,2(r)] ■ exp {-rya^), (4) 

where BK+,i{r) and Bk+,2{i^) are the two degenerate Bloch modes corresponding to the Kj^ 
point, a is the width of the envelope, and J\f determines the beam's power. Surprisingly, 
we find that the diffraction pattern contradicts these expectations by having a three-fold 
symmetry (Fig. El^c)): instead of the circular rings we obtain a triangle pointing to the right. 
Even more surprising is the diffraction pattern obtained when the sign of the nonlinearity 
is reversed: defocusing instead of focusing (Fig. [2t^d)). In the defocusing case, the emerging 
triangular diffraction pattern is almost identical to the focusing case, but is rotated by 
TT. The values of the parameters used in the simulation are D = 35[/im], max{5n} = 
7 ■ 10~^, max{n2|'?/'P} = 5 ■ 10~^, no = 1.5, k = 7.85 ■ 10^[m~^], and the propagation 
distance Z = 8 [cm],. 

We find another interesting phenomenon by repeating the numerical experiment with 
similar input beam that is centered around the other non-equivalent Dirac point K^: the 
orientation of the triangle is reversed again. In other words, for a focusing nonlinearity, 
similar input beams centered around the different Dirac points evolve into triangles that are 
mirror image of each other. We emphasize that the vr rotation of the intensity pattern that 
results from the different Dirac point is exact, as opposed to the vr rotation resulting from 
reversing the sign of the nonlinearity. In the latter case, the orientations of the triangles 
are indeed rotated, but the intensities are not identical, as can be seen from Fig. |2] (c) and 
(d). Therefore, even though similar physical effects are obtained, their physical origin is 
different. 

In all these simulations, we find that the width of the "triangular ring" is sma/Zer than 



4 



the width of the circular ring (arising in the hnear case), for both types (signs) of nonhnear- 
ity. This result is surprising since generally the focusing and the defocusing nonlinearities 
have opposite effects on the width of propagating wave-packet. Namely, in a homogeneous 
medium, the former will always decrease the width of a finite beam, whereas the latter will 
increase it. These opposite tendencies occur also in periodic structures, where the action 
on the beam width is determined by sign of the nonlinearity with respect to the sign of the 
effective mass, e.g., in regions of anomalous diffraction (negative effective mass) a focusing 
nonlinearity always broadens the beam, whereas defocusing narrows it, etc. As stated above, 
the honeycomb lattice is also unique in the sense that the expansion of the evolving trian- 
gular ring-beam is not affected by the sign of nonlinearity; the width of the ring decreases 
for both cases (Fig. M^c) and[2](d)). Moreover, the mean radius of the triangular ring beam, 
defined as 



is the same for the linear and nonlinear diffraction patterns, i.e., the mean radius of the 
evolving beam is not affected by the weak nonlinearity at all. The results of the numerical 
experiment can be summarized as follows. 

1. The nonlinearity transforms conical diffraction into triangular diffraction. 

2. The wave-packet evolves into a triangular ring for both focusing and defocusing non- 
linearity, and nonlinearities of opposite signs result in similar triangles pointing to 
opposite directions. 

3. Two identical wave-packets centered around different Dirac points yet subject to the 
same nonlinearity, evolve to identical triangular rings pointing in opposite directions. 

The remaining part of the paper presents our explanation of these unexpected observa- 
tions that cannot be explained simply by adding a nonlinear term to the Dirac equation. 

In order to analyze the dynamics, we write the wave-function as a linear combination of 
Bloch modes (Bloch decomposition) 



where QmiQi z) is the the amplitude of each Bloch mode, and Bqjri{r) is a Bloch wave defined 
in Eq.(l3]). Substituting the expansion of ip into Eq.([l]) it is clear that the nonlinear term 
mixes four Bloch waves with different lattice momenta, q, i.e., the nonlinear term generates 
spatial four-wave-mixing and therefore the output beam may include of Bloch modes that 
were not initially populated. Hence, even though the input beam is very localized in /c-space 
around the Dirac point, it does not necessarily remain localized during the propagation. In 
fact, previous theoretical and experimental work in nonlinear photonic lattices has clearly 



(5) 
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showed that the Bloch population grows by virtue of four- wave-mixing 28|, and in some 
cases it even leads to spatial super continuum [2^. 

More quantitative analysis of the nonlinear evolution of the distribution of Bloch modes 
reveals the underlying physics of the triangular diffraction pattern. The population of the 
m*'* band is defined as 

P^ = J d'q\{tP{z)\B,,,r.)\' = j (fq\g^{q,z)\\ (7) 

BZ BZ 

and the population imbalance between the first two bands 

A{z) = P,-P, = j d\{\g,{q,z)\'-\g,{q,z)\'). (8) 

BZ 

We calculate the projections of the input and output beams on the Bloch waves of the 
first two bands and find very interesting results: the Bloch distribution |(?m(9;0)p of the 
initial beam has circular symmetry (as expected), with no population imbalance, whereas 
the Bloch distribution of the output beam \gm{qiL)\'^ has three-fold symmetry with a 
significant population imbalance. Fig. |3] presents the evolution of the Bloch distribution for 
the beam of Fig. [2] and the population in each band, for both types of nonlinearity. Clearly, 
focusing nonlinearity increase the population in the first band at the expense of the second 
band, thereby generating a positive imbalance, whereas defocusing nonlinearity generates a 
negative imbalance. While calculating the imbalance, A (2), during the propagation we find 
that it is generated only in the early stages of the propagation (Fig. Hl^a)). Moreover, the 
final imbalance increases with the power of the beam (Fig. ID^b)). 
This analysis leads to three main conclusions: 

1. The nonlinearity gives rise to population transfer to Bloch modes that reside outside 
the Dirac cone, hence, nonlinearity breaks the Dirac dynamics. 

2. The sign of nonlinearity determines the direction of population transfer. 

3. The Bloch distribution of the output beam always has three-fold symmetry in fc-space, 
hence it is clear that the intensity in real space has three-fold symmetry as well. 

These conclusions can be explained by more fundamental arguments, where the most 
crucial step is to understand why the nonlinearity generates a population imbalance between 
the bands. In order to do so, consider the functional 

Hd{z)=HL + 'HNL. (9) 

where 

'Hl = -J iVx^l' - ^Jn{x,yM''^ d\ and (10) 

TiNL^ [ ^n,{\i;\'Yd'x. (11) 
./ -inn 
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Note that Hd is a constant of motion, i.e., dzTid = 0. The evolution equation, Eq.(IT]), is 
derived from T-Ld via the variational principle 

.di^__srLd 

It is instructive to note that I-Ll can be interpreted as the mean propagation constant of the 
wave-packet in the absence of nonlinearity 

nL{z) = {i,{z)\Ho\i^{z)) = [ d'q/3miq)\9m{q,z)\'. (13) 



Since Tid is a constant of motion, we can write 

■Hl{0) + nNLiO) = y-Liz) + UNLiz). (14) 

The initial beam is constructed around the Dirac point with equal population of both bands. 
Since the bands are almost symmetric, from Eq. (fT3l) we find that 'Hl(O) equals the prop- 
agation constant at the Dirac point, Since the spectrum can be shifted by a constant 
without affecting the dynamics, we can set 'Hl(O) = (3d = 0. During the propagation the 
beam experiences significant broadening, and since the total power is conserved the ampli- 
tude of the wave-packet decreases with z. Hence, the nonlinear contribution to Hd becomes 
negligible compared to the linear contribution at large z, i.e., TiNLiz) ^ TiLiz). Therefore 
Eq-dllD simplifies to 

nL{z)c^nNL{0) = [ {m'Yd'x, (15) 

thus the sign of TiLiz) is identical to the sign of n2- By substituting Eq.( !T3|) and using the 
symmetry of the bands, P2iQ) ~ we obtain 

nNLiO) ^ J SqPM) {\gi{q,z)\^ - \g2{q,z)\') . (16) 

BZ 

Since /3i(qr) is positive (first band), the sign of 1-Lnl{S^) is the same as the final population 
imbalance, hence the non-linearity determines the direction of population transfer between 
the bands. From Eqs. f|T5|) and f|T6l) it is clear that positive 77-2 (focusing nonlinearity) yields 
a positive population imbalance, i.e., after a fairly large propagation distance, the first band 
has greater population. For both signs of nonlinearities the resulting Bloch distributions are 
very similar. However, since the group velocity of the bands differs in sign, the orientations 
of the resulting triangles are rotated by tt with respect to each other. We emphasize that, 
other than creating a population imbalance and broadening the Bloch distribution, the 
different signs of nonlinearity have the same physical effect, as opposed to other systems. 
The arguments given above are supported by numerical calculation of 'Hl,'Hnl during the 
propagation, as presented in Fig. [51 Notice that most of the population transfer occurs 
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during the early stages of propagation - just as the population imbalance shown in Fig. 
11(a), and with the same characteristic length. 

Now that we understand that the distribution of Bloch modes varies during propagation, 
we turn to understand the resulting three-fold symmetry. We reexamine the bands around 
the Dirac point by plotting a contour of f3{q) as shown in Fig. [6|, and find that very close to 
the Dirac point, f3{q) has circular symmetry, but further away it has three-fold symmetry. 
Since the new distribution of the Bloch modes follows /3{q), the resulting distribution also 
has three-fold symmetry. This argument explains the Bloch distribution of the output beam 
implying that the intensity in real space also has a three-fold symmetry. Due to ^-reversal 
symmetry, = f3{—q) implying that the contours of equal f3{q) around the point are 
mirror image of the contour lines around the point K_. Therefore, the Bloch distribution 
and hence the real space intensity resulting from the different Dirac cones are exact mirror 
image of each other. 

Next, we try to give more intuitive explanation for the spatial intensity pattern, based 
on the transverse group velocity defined as Vg = ^qP{q), which is in fact the angle of 
propagation. The input beam populates a region in /c-space in which has circular 

symmetry, and all waves propagate with a transverse group velocity of the same magnitude 
in the radial direction. Since the nonlinearity (say, focusing, where we have shown that 
the nonlinear dynamics is determined by the structure of the first band) causes the Bloch 
distribution to expand in /c-space and leave the Dirac cone, some of the modes propagate 
with a greater group velocity than others. After some propagation distance, these modes 
are mapped in real space to the most distant points which are the vertices of the triangle in 
real space. From Fig. [6] it is clear that around the different Dirac points the directions of 
maximal group velocity are opposite. Therefore, a beam centered around with focusing 
nonlinearity evolves into a triangle pointing to the right, whereas the similar input beam 
centered around with focusing nonlinearity evolves into a triangle pointing to the left 
(Fig. int^a)). This explains why similar input beams around the different Dirac points 
subjected to the same nonlinearity evolve into triangles pointing in opposite directions. 

When the same input beam is subjected to nonlinearity of opposite sign, the wave-packet 
expands in fc-space as well. However, for the opposite sign of nonlinearity, the energy is 
transferred to the opposite band (as explained earlier). The contour lines of both bands 
are almost identical, but the group velocity has an opposite sign (Fig. [6]), and therefore the 
corners of the intensity triangle (in real space) are again opposite to each other. Therefore, an 
input beam centered around K+{K-) with focusing non-linearity evolves into a triangular 
ring-beam pointing to the left (right), whereas when the same input beam is subjected to a 
defocusing nonlinearity it evolves into a triangle pointing to the right (left). 

Up to this point the break-down of Dirac dynamics was considered in the context of 
optics. However, we predict the existence of identical phenomena in the context of ultra- 
cold atomic BECs, based on the fact that Eq. ([1]) describes dynamics of an interacting BEC 
in a honeycomb potential. Such potentials can be constructed optically, however, they were 
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considered mainly in the context of fermionic gases (e.g., see Refs. |5|-|8|), probably because 
electrons in real graphene are fermions, for which the ground state is a filled Fermi sea up the 
Fermi level (which can be shifted from Dirac point depending on doping), and because the 
ground state of a Bose gas in a potential with graphene band structure is not a condensate 
occupying Bloch modes in the vicinity of the Dirac point. Thus, in order to obtain such a 
state one should do it dynamically. Here we propose a scheme which is not trivial, but still 
realistic experimentally: Suppose that one prepares a finite size condensate in a harmonic, 
pancake shaped, effectively two-dimensional trap 30|. At some point in time, the in-plane 
confinement is turned off (the confinement keeping the 2D structure is still present), the 
pancake EEC is subsequently given a momentum kick with three plane waves with angles 
27r/3 with respect to each other, which are directed into the three equivalent corners of the 
BZ (say -R'+), and that the honeycomb potential is ramped up (the lattice spacing of the 
honeycomb potential should be much smaller than the size of the cloud); the initial EEC can 
be noninteracting, and the nonlinearity can be turned on via the Feshbach mechanism (say 
when the honeycomb potential is ramped up). Such an experimental sequence of events is not 
an easy task but it is within today's experimental capabilities. The sign of the nonlinearity 
and the strength of the interactions can be tuned for some EECs by a Feshbach resonance 
mechanism. 

In conclusion, honeycomb lattices possess intersection points between the first two bands, 
and Eloch modes that reside close to these points give rise to extremely unique dynamics: 
dynamics of massless spin half particles described by the Dirac equation. We have shown 
that, in the presence of nonlinearity, this unique dynamics is significantly altered, and the 
Dirac equation is no longer suitable for describing the (propagation) evolution. In addition, 
we studied the nonlinear evolution of wave-packets comprised of such Eloch modes, which 
under linear conditions yield conical diffraction. We have found that both signs of nonlin- 
earity transform the circular rings into triangular rings, resulting in triangular diffraction 
pattern. This new type of diffraction cannot be obtained simply by adding the nonlinear 
term into the effective Dirac equation, since the Eloch modes distribution quickly expands 
beyond the range of the Dirac cone, where the effective evolution equation governing the 
dynamics cannot be approximated by the Dirac equation, nor by its nonlinear extension. 
Hence we expect that the intriguing phenomena obtained in honeycomb lattices such as 



Klein tunneling, Zitterbewegung, anti-localization, zero modes edge-states [3l|, |32| and more 
would be significantly altered by the presence of nonlinearity. Likewise, the intriguing fea- 
tures related to Dirac points in photonic crystals j9-12| would be affected by nonlinearities 
of the same type studied here. 
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FIG. 1: (a) Perfect honeycomb lattice that has two sites in a unit cell, (b) The first Brillouin zone 
with the high symmetry points, (c) The first two bands, and (d) The vicinity of the intersection 
(Dirac) point. 
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FIG. 2: (a) Input beam from the vicinity of K^. (b) The output of the linear propagation, (c) the 
output of the nonlinear propagation with focusing nonlinearity, and (d) the output of the nonlinear 
propagation with defocusing nonlinearity. 
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FIG. 3: The bloch mode decomposition of the input wave-packet from the vicinity of the Dirac 
point, (a)-(d) The Bloch distribution of the first band and (e)-(h) the Bloch distribution 

of the second band |(72(q)P- In each figure we write the corresponding population. Clearly, focus- 
ing nonlinearity transfers the energy to the first band, and defocusing non-linearity transfers the 
population to the second band. 
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FIG. 4: (a) The population in the first two bands during the propagation in the presence of 
defocusing nonlinearity, where max{|n2| IV'P} = 4 • 10~^. (b) The final population imbalance as a 
function of the strength of the nonlinear refractive index. 
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FIG. 5: (a) Calculation of T-Lci{z) and in the presence focusing nonlinearity. (b) The nonlinear 

contribution to 7ici{z), that decreases rapidly during the propagation. 
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FIG. 6: Contour plot of the first two bands. The 
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Dirac point are denoted by K±. The black arrows 
is maximal. 
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